Figure 4 and supplementary figure 4

Author

Niccolò Arecco

Last code execution: 2024 February 07, Wednesday @ 19:02:17.

1 Intro

We performed histone PTMs DIA mass spectrometry and I analysed the raw files with EpiProfile v2.1 (Yuan et al. 2018). Here I parse the EpiProfile output and analyse the PTMs values with DEP (Zhang et al. 2018).

Acknowledgements: I’d like to thank Zuo-Fei Yuan and Simone Sidoli from the Ben Garcia lab for their help with data analysis and sample preparation respectively. Also a big thank you to Cristina Chiva from the CRG/UPF proteomics facility.

2 Set Up

To fetch files on the CRG cluster I mount the server on my local computer with the following command in the terminal.

Code
sshfs narecco@ant-login.linux.crg.es:/users/mirimia ~/mnt

2.1 Packages

Load packages required for the analysis and suppress any message. Check the Section 6 section at the end for package versions.

Code
library(dplyr, warn.conflicts = F, quietly = T)
library(readr)
library(readxl)
library(pheatmap)
library(DEP)
library(Cairo)
library(DT)

Load my R package niar to use my custom made function to parse EpiProfile output.

Code
if ( ('niar' %in% .packages(all.available = TRUE)) == TRUE ) {
  library('niar')
} else if ( ('niar' %in% .packages(all.available = TRUE)) == FALSE ){
  message("The package niar is not installed so I'm trying to load it from ",
          "the local repository of the package")
  if ( dir.exists('~/mnt/narecco/software/R/niar' ) )  {
    devtools::load_all(path = '~/mnt/narecco/software/R/niar')
  } else {
    stop("Can't find the local repo of the niar package! ",
         "You must install it with:\n",
         "devtools::install_github('Ni-Ar/niar') ")
  }
} else{
  stop("Can't understand if 'niar' package was installed beforehand")
}
The package niar is not installed so I'm trying to load it from the local repository of the package
ℹ Loading niar

2.2 Functions

Define a plot style.

Code
theme_classic(base_size = 6.3, base_family = "Arial") +
  theme(plot.title = element_text(size = 8, vjust = -1, hjust = 0),
        plot.background = element_blank(),
        strip.background = element_blank(),
        panel.background = element_blank(),
        axis.line = element_line(linewidth = 0.1),
        axis.title.x = element_blank(),
        axis.title.y = element_text(size = 5),
        axis.text = element_text(colour = "black"),
        axis.text.x = element_text(margin = margin(t = 0, unit = "mm")),
        axis.ticks.x = element_blank(),
        axis.ticks.y = element_line(linewidth = 0.1),
        axis.ticks.length.y = unit(1, "mm"),
        legend.margin = margin(l = 0, t = 0, unit = "mm"),
        legend.box.margin = margin(l = 0, unit = "mm"),
        legend.box.background = element_blank(),
        legend.box.spacing = unit(0.5, "mm"),
        # legend.title = element_blank(),
        legend.text = element_text(margin = margin(l = -10, unit = "mm")),
        legend.text.align = 0,
        legend.key.height = unit(2, "mm"),
        legend.key.width = unit(7.5, "mm"),
        legend.key = element_blank(),
        legend.background = element_blank(),
        legend.spacing.x = unit(3, 'mm')
        ) -> alluvial_thm 

I modify the geom_stratum() function just to add a very small line width the the black borders of the columns. With linewidth = .05, I don’t have to change them manually later in Illustrator. This will hopefully fixed in future

Small function to parse the PTM abundance as 10 to the power scientific notion.

Code
scientific_10 <- function(x) {
  parse(text=gsub("e\\+", " %*%10^", scales::scientific_format()(x)))
}

2.3 Directories & File Paths

Here I organise all the file and directories paths I need to run the analysis and define where to save the processed tables and figures.

Code
oneDrive_Dir <- file.path("~/OneDrive - CRG - Centre de Regulacio Genomica/Suz12_AS_project")  
code_dir_fig <- file.path(oneDrive_Dir, "_Code/Fig4")
tbl_dir_fig4 <- file.path(code_dir_fig, "tables")
pdf_dir_fig4 <- file.path(code_dir_fig, "pdfs")

if (!dir.exists(pdf_dir_fig4)) { dir.create(pdf_dir_fig4, recursive = T) }
if (!dir.exists(tbl_dir_fig4)) { dir.create(tbl_dir_fig4, recursive = T) }

Path to medata data and EpiProfile output files.

Code
metadata_path <- file.path(tbl_dir_fig4, 'Histone_PTMs_metadata.xlsx')

WT_Dex4_histone_ratio_path <- list.files(tbl_dir_fig4,
                                         pattern = "FullSet1_Run2_DIA_histone_ratios_Suz12dExon4.xls",
                                         full.names = T)

Resce_KO_histone_ratio_path <- list.files(tbl_dir_fig4,
                                          pattern = "FullSet2_Run1_DIA_histone_ratios_Suz12KO.xls",
                                          full.names = T)

stopifnot(file.exists(metadata_path))
stopifnot(file.exists(WT_Dex4_histone_ratio_path))
stopifnot(file.exists(Resce_KO_histone_ratio_path))

SUZ12 WT proteoforms quantification data.

Code
protein_psi_path <- file.path(tbl_dir_fig4, 'Protein_PSI_3reps_LS.xlsx')

3 Main Figure Panels

Import metadata.

Code
mtdt <- read_excel(metadata_path) |> 
  setNames(c("Sample", "Short_Name", "Clone_Name", "Condition", "Replicate"))

mtdt$Condition <- factor(mtdt$Condition, 
                         levels = c("WT", "CSexon4", "Dexon4", "KO", "KO_L", "KO_S") )

3.1 H3 K27-K40 peptide modification ratios

Import DIA data quantification using automatic EpiProfile v2.1 quantifications.

Code
read_EpiProfile_histone_ratio(WT_Dex4_histone_ratio_path) |>
  left_join(y = mtdt, by = "Sample") |>
  mutate(Facility_ID = Sample) |>
  mutate(Sample = Short_Name) |>
  tidy_hPTMs(split_double_PTMs = T) -> dia_auto_set1

read_EpiProfile_histone_ratio(Resce_KO_histone_ratio_path) |>
  # fix an extra _ in the file name
  mutate(Sample = gsub(pattern = '2022MH003_NIAR_005_02__30pto_DIA', 
                       replacement = '2022MH003_NIAR_005_02_30pto_DIA', 
                       x = Sample) ) |>
  left_join(y = mtdt, by = "Sample") |>
  mutate(Facility_ID = Sample) |>
  mutate(Sample = Short_Name) |>
  tidy_hPTMs(split_double_PTMs = T) -> dia_auto_set2

data_dia_auto <- rbind(dia_auto_set1, dia_auto_set2)

Tidy up the information for histone peptide H3 K27- K36. Here I split the H3 K27 and K36 signal coming form the same peptide. In this way I can pile up the all methyl-forms.

Code
data_dia_auto |>
  subset(Peptide_Start >= 27 & Peptide_End <= 40) |>
  mutate(K27 = case_when( is.na(Sub_PTM1) & grepl(pattern = "^K27", x = Modification) ~ Modification,
                          grepl(pattern = "unmod", x = Modification) ~ "unmod" ) ) |>
  mutate(K36 = case_when( is.na(Sub_PTM2) & grepl(pattern = "^K36", x = Modification) ~ Modification, 
                          grepl(pattern = "unmod", x = Modification) ~ "unmod" ) ) |>
  relocate(K27, K36, .after = Sub_PTM2) |>
  unite(col = "K27", c("Sub_PTM1", "K27"), sep = "", na.rm = T, remove = T ) |>
  unite(col = "K36", c("Sub_PTM2", "K36"), sep = "", na.rm = T, remove = T ) |>
  mutate(K27 = case_when(grepl(pattern = "^$", x = K27) ~ "unmod",
                         grepl(pattern = "[aA-zZ]", x = K27) ~ K27 )) |>
  mutate(K36 = case_when(grepl(pattern = "^$", x = K36) ~ "unmod",
                         grepl(pattern = "[aA-zZ]", x = K36) ~ K36 )) |>
  # When one mark is not found I label it unmodified ('unmod') but it could also
  # be labelled as "only the other mark".
  # mutate(K27 = case_when(grepl(pattern = "^$", x = K27) ~ "Only\nK36me",
  #                        grepl(pattern = "[aA-zZ]", x = K27) ~ K27 )) |>
  # mutate(K36 = case_when(grepl(pattern = "^$", x = K36) ~ "Only\nK27ac-me",
  #                        grepl(pattern = "[aA-zZ]", x = K36) ~ K36 )) |>
  relocate(K27, K36, .after = PTM) |>
  mutate(Histone = factor(Histone, levels = c("H3", "H3.3") ) ) -> data_K27_36

Overview of H3 and H3.3 PTMs in lysines 27 and 36

Code
datatable(data_K27_36, rownames = FALSE, filter = 'top', 
          options = list(pageLength = 7, autoWidth = TRUE) )

I can also deconvolve the K27 and K36 signals from the same peptide and pile-up each signal.

Code
data_K27_36 |>
  group_by(Condition, K27, Histone) |>
  mutate(Average_Ratio = mean(Ratio)) |>
  mutate(K27 = factor(K27, levels = c("K27ac", "K27me3", "K27me2", "K27me1", "unmod") ) )|>
  select(Histone, Condition, Average_Ratio, K27, Short_Name, Ratio) |>
  mutate(Ratio_Sum = sum(Average_Ratio) / 3) |>
  mutate(Percent = round(Average_Ratio*100, 2) ) |>
  select(Condition, Histone, K27, Average_Ratio, Ratio_Sum, Percent) |>
  unique() |>
  mutate(across(.cols = where(is.numeric), \(x) round(x, 4) ) ) |>
  arrange(Condition, Histone, K27) |> ungroup() -> data_K27

Same as above but for K36.

Code
data_K27_36 |>
  group_by(Condition, K36, Histone) |>
  mutate(Average_Ratio = mean(Ratio)) |>
  mutate(K36 = factor(K36, levels = c("K36me3", "K36me2", "K36me1", "unmod") ) ) |> 
  select(Histone, PTM, Condition, Average_Ratio, K36, Short_Name, Ratio) |>
  mutate(Ratio_Sum = sum(Average_Ratio) / 3) |>
  mutate(Percent = round(Average_Ratio*100, 2) ) |>
  select(Condition, Histone, K36, Average_Ratio, Ratio_Sum, Percent) |>
  unique() |>
  mutate(across(.cols = where(is.numeric), \(x) round(x, 4) ) ) |>
  arrange(Condition, Histone, K36) |> ungroup() -> data_K36

For specific histone PTMs one can explore the results in the following tables.

Overview of H3 and H3.3 K27 acetylation and methylation quantifications

Code
datatable(data_K27, rownames = FALSE, filter = 'top', 
          options = list(pageLength = 5, autoWidth = TRUE) )

Overview of H3 and H3.3 K36 methylation quantifications

Code
datatable(data_K36, rownames = FALSE, filter = 'top', 
          options = list(pageLength = 5, autoWidth = TRUE) )

3.1.1 Global H3K27 ratios

Plot figure 4B panel with the pileup of figure H3K27 methylation and acetylation.

Code
data_K27_36 |>
  group_by(Condition, K27, Histone) |>
  mutate(Average_Ratio = mean(Ratio)) |>
  summarise(Ratio_Sum = sum(Average_Ratio) / 3) |> # 3 replicates
  subset(Histone %in% "H3") |>
  mutate(K27 = factor(K27, levels = c("K27ac", "K27me3", "K27me2", "K27me1", "unmod") ) ) |>
  ggplot() +
  aes(x = Condition, stratum = K27, alluvium = K27, 
      y = Ratio_Sum, fill = K27) +
  geom_alluvium(alpha = 0.5, width = 0.4) +
  geom_stratum2(alpha = 1, width = 0.6, colour = "black") +
  geom_fit_text(aes(label = K27 ),
                stat = "stratum", width = 0.6, min.size = 2,
                padding.x = grid::unit(0.2, "mm"),
                padding.y = grid::unit(0.2, "mm"),
                size = 7, family = "Arial", show.legend = F) +
  scale_fill_manual(values = c("K27ac" = "#feac81", "K27me1" ="#ced1af", 
                               "K27me2" = "#748f46", "K27me3" = "#47632a", 
                               "unmod" = "gray" ),
                    name = 'H3') +
  scale_x_discrete(expand = expansion(mult = 0.05, add = c(0, 0.1)),
                   labels = c('WT', 'CSex4', '∆ex4', 'KO', 
                              'KO+L', 'KO+S')) +
  scale_y_continuous(expand = expansion(mult = 0, add = 0) ) +
  guides(fill = guide_legend(byrow = TRUE) ) + # required for the legend item spacing 
  labs(y = "Average Ratio") +
  alluvial_thm -> p_Alluvial_K27

p_Alluvial_K27 

Save to pdf.

Code
ggsave(filename = "Fig4B_H3_only_K27_Mrgd_Alluvial.pdf", path = pdf_dir_fig4, 
       plot = p_Alluvial_K27, device = cairo_pdf, units = "cm",
       width = 6.5, height = 4.85)

4 Differential abundance analysis of histone PTMs with DEP

Use the area under the peak as counts to be analysed using DEP framework.

Code
rbind( read_EpiProfile_histone_ratio(WT_Dex4_histone_ratio_path),
       read_EpiProfile_histone_ratio(Resce_KO_histone_ratio_path) ) |>
  # fix an extra _ in the file name like before
  mutate(Sample = gsub(pattern = '2022MH003_NIAR_005_02__30pto_DIA', 
                      replacement = '2022MH003_NIAR_005_02_30pto_DIA', 
                      x = Sample) ) |>
  left_join(y = mtdt, by = "Sample") |>
  mutate(Facility_ID = Sample) |>
  mutate(Sample = Short_Name) -> data_dia_auto_not_tidy

One small caveat in this analysis is that the areas under the peaks are very large numbers or at least larger than 2147483647 which is the highest number R can store as integer. To convert the areas to integers from numeric type one could use the package bit64 that stores integers in 64 bits instead of the default 32 bit. My solution was however simpler, I divided the all values by 1000 to scale down by 3 order of magnitude. I then convert the dataframe to matrix.

I tried to add a pseudocounts to the data but it actually doesn’t help with the analysis and is better to input the areas as they are. The function make_unique is the first step of the analysis of DEP and creates a unique data frame.

Code
data_dia_auto_not_tidy |>
  tidy_hPTMs() |>
  select(First_Col, PTM, Sample, Area) |>
  # reduce all areas by a 1000 to avoid having to coerc to integer big numbers
  mutate(Area = Area / 1e4) |>
  pivot_wider(id_cols = c(First_Col, PTM), names_from = Sample, values_from = Area) |>
  # Coerce the areas to integers
  mutate(across(.cols = ends_with(c("_1", "_2", "_3")), .fns = as.integer )) |>
  make_unique(names = "First_Col", ids = "PTM") -> data_unique

The data frame data_unique contains 22 columns (of which 18 are samples) and 185 rows corresponding to the identified histone PTMs.

Create a metadata compatible for DEP set the KO as the reference sample as I’m interesting in the rescue efficiency relative to the KO.

Code
# Set metadata colnames to be as `make_se` wants them to be 
colnames(mtdt) <- c("Sample", "label", "Clone_Name", "condition", "replicate")

mutate(mtdt, condition = factor(condition, levels = c("WT", "CSexon4", "Dexon4",
                                                      "KO", "KO_L", "KO_S"))) |>
  print(n = 3) -> mtdt_dep
# A tibble: 18 × 5
  Sample                          label Clone_Name condition replicate
  <chr>                           <chr> <chr>      <fct>         <dbl>
1 2022MZ006_IVMO_001_02_25pto_DIA WT_1  ZWT        WT                1
2 2022MZ006_IVMO_002_02_25pto_DIA WT_2  ZA1        WT                2
3 2022MZ006_IVMO_003_02_25pto_DIA WT_3  ZA2        WT                3
# ℹ 15 more rows

Find the columns that contain the samples areas in the data_unique dataframe and create a SummarizedExperiment object to be used in for the DEP analysis. The first step is to do a log2 transformation of the area counts.

Code
Samples_area_columns <- grep(pattern = paste(mtdt_dep$label, collapse = "|"), colnames(data_unique))
data_se <- make_se(data_unique, columns = Samples_area_columns, expdesign = mtdt_dep)

Now the data_se contains the modifications in the rows (assuming they are proteins) and in the column the area.

Code
data_filt <- filter_missval(data_se, thr = 0)

Check PTMs frequency across samples.

Code
plot_frequency(data_filt)

Overall most modifications are found in all samples.

Code
data_norm <- normalize_vsn(data_filt)

Verify the fit.

Code
meanSdPlot(x = data_norm)

Seems okay.

Check normalisation effect on data.

Code
plot_normalization(data_se, data_filt, data_norm) +
  scale_fill_manual(values = c("WT" = "royalblue3", "Dexon4" = "goldenrod",
                               "CSexon4" = "#B0461C", "KO" =  "gray45", 
                               "KO_L" = "mediumpurple3", "KO_S" = "darkorange1")) +
  theme(axis.text = element_text(color = 'black'),
        plot.background = element_blank(), 
        panel.background = element_blank())

Check missing values.

Code
plot_missval(data_norm)

I don’t see obvious missing values trends, besides only that replicate #3 of CSexon has less PTMs detected. I don’t think it’s necessary to impute missing values.

4.1 Test for significance against WT control

Test every sample versus WT control

Code
data_diff <- test_diff(data_norm, type = "control", control = "WT")

Significance thresholds:

  • P-value <= 0.05

  • |Log2 Fold Change| >= 1.5

Code
dep <- add_rejections(data_diff, alpha = 0.05, lfc = 1.5 )

Check sample clustering a Gower’s distance matrix.s

Code
plot_dist(dep, significant = TRUE, pal = "PuOr") 

This is basically showing the 2 batches of the 2 MS acquisition experiments.

Check sample clustering using the same distance.

Code
plot_heatmap(dep, type = "centered", kmeans = F, clustering_distance = "gower",
             col_limit = 4, show_row_names = TRUE,
             indicate = c("condition", "replicate"))

Get results to make a volcano plot.

Code
data_results <- get_results(dep)
data_results |>
  as_tibble() |>
  select(-significant) |>
  pivot_longer(cols = !c("name", "ID"),
               names_to = c("contrast", "Feature"),
               names_sep = c("_(vs_WT_|centered)"),
               values_to = "Value"
               ) |>
  mutate(Feature = ifelse(test = Feature == "", yes = "centered", no = Feature)) |>
  pivot_wider(names_from = "Feature", values_from = "Value") |>
  dplyr::select(-ID) -> res_all_tidy

colnames(res_all_tidy)[colnames(res_all_tidy) == "name"] <- 'First_Col'

Extract WT sample area of each PTM, to be used for the volcano plot.

Code
data_dia_auto |>
  subset(Condition == "WT") |>
  group_by(First_Col) |>
  mutate(WT_Mean_Area = mean(Area) + 1 ) |>
  relocate(WT_Mean_Area, .after = PTM) |>
  select(c(First_Col, WT_Mean_Area)) |>
  unique() -> WT_area

Get the histone PTMs info columns.

Code
histone_ptms_anno <- unique(data_dia_auto[, 1:10])

Add all the info the the tidy results dataframe.

Code
res_all <- left_join(x = res_all_tidy, y = histone_ptms_anno, by = join_by(First_Col))
res_all <- left_join(res_all, WT_area, by = join_by(First_Col) )

res_all |>
  mutate(Histone = factor(Histone, levels = sort(unique(res_all$Histone)) ) ) |>
  arrange(desc(Histone), Peptide_Start, Peptide_End) |>
  ungroup() -> res_all

Export all tidy results for supplementary table 3 as tab delimited, or as csv with UTF-8 Byte order mark which indicates to Excel the csv is UTF-8 encoded.

Code
write_delim(x = res_all, file = file.path(tbl_dir_fig4, 'DEP_analysis_histone_MS.tab'), 
            delim = '\t', col_names = T, append = F, quote = 'none', na = "NA",
            progress = F, escape = 'none')

write_excel_csv(x = res_all, 
                file = file.path(tbl_dir_fig4, 'DEP_analysis_histone_MS.csv'),
                 delim = ',', col_names = T, append = F, quote = 'none', 
                na = "NA", progress = F, escape = 'none')

Interactively explore the table

Code
res_all |>
  mutate(across(.cols = where(is.numeric), \(x) round(x, 4) ) ) |>
  datatable(rownames = FALSE, filter = 'top', 
            options = list(pageLength = 3, autoWidth = TRUE) )

4.2 PTMs centred abundances

Plot centred abundance

Code
res_all |>
  subset(Peptide_Start >= 27 & Peptide_End <= 40) |>
  subset(!is.na(centered)) |>
  ggplot(aes(x = centered, y = PTM, fill = contrast)) +
  facet_grid(rows = vars(Modification), scales = 'free' ) +
  geom_vline(xintercept = 0, linewidth = 0.4) + 
  geom_col(width = 0.75, 
           position = position_dodge(width = 0.75, preserve = 'single')) +
  geom_point(aes(size = WT_Mean_Area),
             shape = 21, position = position_dodge(width = 0.75),
             stroke = 0.2) +
  labs(x = "Centered PTM abundance") +
  scale_x_continuous(limits = c(-5, 5), oob = scales::oob_squish ) +
  scale_size_area(breaks = c(1e8, 1e9, 1e10, 1e11, 1e12),
                  labels = scientific_10, name = "WT PTM\nabundance") +
  scale_fill_manual(values = c("WT" = "royalblue3", "Dexon4" = "goldenrod",
                               "CSexon4" = "#B0461C", "KO" =  "gray45", 
                               "KO_L" = "mediumpurple3", "KO_S" = "darkorange1"),
                    name = "Sample") +
  theme_classic() +
  theme(strip.text = element_blank(),
        strip.background = element_blank(),
        axis.text = element_text(colour = "black"),
        axis.title.y = element_blank(),
        legend.key.size = unit(1, 'mm'),
        panel.grid.major.y = element_line()) -> p_centered_K27_K36

p_centered_K27_K36

4.3 Volcano plot

Show differential histone PTMs as volcano plots.

4.3.1 ∆ex4 vs WT

Prepare data for plot.

Code
res_Dex4 <- subset(res_all, contrast == "Dexon4" & !is.na(p.val))

res_Dex4 <- res_Dex4 |>
  mutate(Direction = case_when(p.adj <= 0.05 & ratio > 0 ~ 'UP',
                               p.adj <= 0.05 & ratio < 0 ~ 'DOWN',
                               p.adj > 0.05 ~ 'None') ) 

Plot figure 4C.

Code
ggplot(res_Dex4 ) +
  aes(x = ratio, y = -log10(p.val), fill = Direction, size = WT_Mean_Area ) +
  geom_point(shape = 21, stroke = 0.2) +
    annotate(geom = "label", x = 5.5,  y = 0.10, label = "∆ex4", 
           colour = "black", fill = "goldenrod", size = 2,
           label.padding = grid::unit(0.5, "mm"), 
           label.r = unit(0.25, "mm"), family = "Arial",
           label.size = grid::unit(0.125, "mm") ) +
  annotate(geom = "label", x = -5.5, y = 0.10, label = "WT", 
           colour = "white",  fill = "royalblue3", size = 2, 
           label.padding = grid::unit(0.5, "mm"), 
           label.r = unit(0.25, "mm"), family = "Arial",
           label.size = grid::unit(0.125, "mm") ) +
  # points on the right side
  geom_label_repel(data = subset(res_Dex4, p.adj < 0.09 & ratio > 0 ), 
                   aes(label = PTM), fill = 'white',
                   seed = 16, show.legend = F, segment.curvature = -1e-20, 
                   family = "Arial", size = 1.85, nudge_x = -0.25,
                   segment.color = 'black',verbose = F, 
                   box.padding = grid::unit(1, "mm"),
                   point.padding = grid::unit(0.55, "mm"),
                   label.padding = grid::unit(0.5, "mm"), 
                   label.size = 0.120, 
                   max.overlaps = 20)  +
  # points on the left side
  geom_label_repel(data = subset(res_Dex4,  p.adj < 0.1 & ratio < 0 ), 
                 aes(label = PTM), fill = 'white',
                 seed = 16, show.legend = F, segment.curvature = -1e-20, 
                 family = "Arial", size = 1.85, nudge_x = 0.25,
                 segment.color = 'black',verbose = F, 
                 box.padding = grid::unit(1, "mm"),
                 point.padding = grid::unit(0.55, "mm"),
                 label.padding = grid::unit(0.5, "mm"), 
                 label.size = 0.120, 
                 max.overlaps = 20)  +
  labs(x = expression(log[2] ~ "Fold Change"), y = expression(-log[10] ~ "P-Value")) +
  scale_fill_manual(values = c('dodgerblue', "gray84", "firebrick3"), guide = 'none') +
  guides(size = guide_legend(override.aes = list(fill = "white"))) +
  scale_x_continuous(expand = expansion(add = 0.02, mult = 0), limits = c(-6.4, 6.4), n.breaks = 7) +
  scale_y_continuous(expand = expansion(add = c(0, 0.25), mult = 0), limits = c(0, NA)) +
  scale_size_area(breaks = c(1e8, 1e9, 1e10, 1e11, 1e12),
                  labels = scientific_10, name = "WT PTM\nabundance") +
  theme_classic(base_size = 6, base_family = "Arial") +
  theme(panel.background = element_blank(),
        panel.border = element_blank(),
        panel.grid.major = element_line(linewidth = 0.2),
        legend.position = c(0.91, 0.87), 
        legend.key = element_blank(),
        legend.key.size = unit(1.0, units = 'mm'),
        legend.text = element_text(size = 5, margin = margin(l = -1, r = -1, unit = 'mm')),
        legend.direction = "vertical",
        legend.background = element_blank(),
        legend.box.background = element_blank(),
        legend.box.spacing = unit(-3, 'mm'),
        axis.text = element_text(colour = 'black'),
        axis.ticks.length = unit(1, units = 'mm'),
        axis.ticks = element_line(colour = 'black'),
        axis.title = element_text(size = 5),
        axis.line = element_line(linewidth = 0.2),
        plot.background = element_blank() ) -> p_Dex4

p_Dex4

Save to pdf.

Code
ggsave(path = pdf_dir_fig4, filename = "Fig4C_Dex4_vs_WT_Volcano_w_labels.pdf",
       plot = p_Dex4, device = cairo_pdf, units = "cm",
       height = 5.4, width = 7.0)

4.3.2 CSex4 vs WT

Prepare data of CSex4 vs WT

Code
res_CSex4 <- subset(res_all, contrast == "CSexon4" & !is.na(p.val))

res_CSex4 <- res_CSex4 |>
  mutate(Direction = case_when(p.adj <= 0.05 & ratio > 0 ~ 'UP',
                              p.adj <= 0.05 & ratio < 0 ~ 'DOWN',
                              p.adj > 0.05 ~ 'None') ) 

Note: nothing is differentially expressed in this contrast

Code
table(res_CSex4$Direction)

None 
 159 

4.3.3 Volcano plot KO vs WT

Code
res_KO <- subset(res_all, contrast == "KO" & !is.na(p.val))

res_KO <- res_KO |>
  mutate(Direction = case_when(p.adj <= 0.05 & ratio > 0 ~ 'UP',
                              p.adj <= 0.05 & ratio < 0 ~ 'DOWN',
                              p.adj > 0.05 ~ 'None') ) 

Plot KO vs WT

Code
ggplot(res_KO ) +
  aes(x = ratio, y = -log10(p.val), fill = Direction, size = WT_Mean_Area ) +
  geom_point(shape = 21, stroke = 0.2) +
    annotate(geom = "label", x = 8.5,  y = 0.10, label = "KO", 
           colour = "white", fill = "gray16", size = 2,
           label.padding = grid::unit(0.5, "mm"), 
           label.r = unit(0.25, "mm"), family = "Arial",
           label.size = grid::unit(0.125, "mm") ) +
  annotate(geom = "label", x = -12.5, y = 0.10, label = "WT", 
           colour = "white",  fill = "royalblue3", size = 2, 
           label.padding = grid::unit(0.5, "mm"), 
           label.r = unit(0.25, "mm"), family = "Arial",
           label.size = grid::unit(0.125, "mm") ) +
  # points on the right side
  geom_label_repel(data = subset(res_KO, p.adj <= 0.15 & ratio > 0 ), 
                   aes(label = PTM), fill = 'white',
                   seed = 16, show.legend = F, segment.curvature = -1e-20, 
                   family = "Arial", size = 1.85, nudge_x = -0.25,
                   segment.color = 'black',verbose = F, 
                   box.padding = grid::unit(1, "mm"),
                   point.padding = grid::unit(0.55, "mm"),
                   label.padding = grid::unit(0.5, "mm"), 
                   label.size = 0.120, 
                   max.overlaps = 20)  +
  # points on the left side
  geom_label_repel(data = subset(res_KO, p.adj <= 0.15 & ratio < 0 ), 
                 aes(label = PTM), fill = 'white',
                 seed = 16, show.legend = F, segment.curvature = -1e-20, 
                 family = "Arial", size = 1.85, nudge_x = 0.25,
                 segment.color = 'black',verbose = F, 
                 box.padding = grid::unit(1, "mm"),
                 point.padding = grid::unit(0.55, "mm"),
                 label.padding = grid::unit(0.5, "mm"), 
                 label.size = 0.120, 
                 max.overlaps = 20)  +
  labs(x = expression(log[2] ~ "Fold Change"), y = expression(-log[10] ~ "P-Value")) +
  scale_fill_manual(values = c('dodgerblue', "gray84", "firebrick3"), guide = 'none') +
  guides(size = guide_legend(override.aes = list(fill = "white"))) +
  scale_x_continuous(expand = expansion(add = 0.02, mult = 0), limits = c(-15, 10), n.breaks = 7) +
  scale_y_continuous(expand = expansion(add = c(0, 0.25), mult = 0), limits = c(0, NA)) +
  scale_size_area(breaks = c(1e8, 1e9, 1e10, 1e11, 1e12),
                  labels = scientific_10, name = "WT PTM\nabundance") +
  theme_classic(base_size = 6, base_family = "Arial") +
  theme(panel.background = element_blank(),
        panel.border = element_blank(),
        panel.grid.major = element_line(linewidth = 0.2),
        legend.position = c(0.91, 0.67), 
        legend.key = element_blank(),
        legend.key.size = unit(1.0, units = 'mm'),
        legend.text = element_text(size = 5, margin = margin(l = -1, r = -1, unit = 'mm')),
        legend.direction = "vertical",
        legend.background = element_blank(),
        legend.box.background = element_blank(),
        legend.box.spacing = unit(-3, 'mm'),
        axis.text = element_text(colour = 'black'),
        axis.ticks.length = unit(1, units = 'mm'),
        axis.ticks = element_line(colour = 'black'),
        axis.title = element_text(size = 5),
        axis.line = element_line(linewidth = 0.2),
        plot.background = element_blank() ) -> p_KO

p_KO

Note: these WT and KO samples were processed at different moment, so there could be some batch effect difference that limits the statistical power.

Save to pdf.

Code
ggsave(path = pdf_dir_fig4, filename = "FigExtra_KO_vs_WT_Volcano_w_labels.pdf",
       plot = p_KO, device = cairo_pdf, units = "cm",
       height = 12, width = 10)

4.4 Test for significance vs the KO+L

Test against Long Suz12 isoform KO rescue.

Code
data_diff_long <- test_diff(data_norm, type = "control", control = "KO_L")

Use same significance thresholds as before.

Code
dep_rescues <- add_rejections(data_diff_long, alpha = 0.05, lfc = 1.5 )

Get results to make a volcano plot of the long vs the short

Code
data_results <- get_results(dep_rescues)
data_results |>
  as_tibble() |>
  select(-significant) |>
  pivot_longer(cols = !c("name", "ID"),
               names_to = c("contrast", "Feature"),
               names_sep = c("_(vs_KO_L_|centered)"),
               values_to = "Value"
               ) |>
  mutate(Feature = ifelse(test = Feature == "", yes = "centered", no = Feature)) |>
  pivot_wider(names_from = "Feature", values_from = "Value") |>
  dplyr::select(-ID) -> res_all_tidy

colnames(res_all_tidy)[colnames(res_all_tidy) == "name"] <- 'First_Col'

Add all the info the the tidy results dataframe.

Code
res_all <- left_join(x = res_all_tidy, y = histone_ptms_anno, by = join_by(First_Col))
res_all <- left_join(res_all, WT_area, by = join_by(First_Col) )

res_all |>
  mutate(Histone = factor(Histone, levels = sort(unique(res_all$Histone)) ) ) |>
  arrange(desc(Histone), Peptide_Start, Peptide_End) |>
  ungroup() -> res_all

Prepare data for plot.

Code
res_KOrS <- subset(res_all, contrast == "KO_S" & !is.na(p.val))

res_KOrS <- res_KOrS |>
  mutate(Direction = case_when(p.adj <= 0.05 & ratio > 0 ~ 'UP',
                               p.adj <= 0.05 & ratio < 0 ~ 'DOWN',
                               p.adj > 0.05 ~ 'None') ) 

4.4.1 Volcano plot rescue SUZ12-L vs SUZ12-S

Plot rescues.

Code
ggplot(res_KOrS) +
  aes(x = ratio, y = -log10(p.val), fill = Direction, size = WT_Mean_Area ) +
  geom_point(shape = 21, stroke = 0.2) +
    annotate(geom = "label", x = 7.5,  y = 0.10, label = "KO+S", 
           colour = "black", fill = "darkorange1", size = 2,
           label.padding = grid::unit(0.5, "mm"), 
           label.r = unit(0.25, "mm"), family = "Arial",
           label.size = grid::unit(0.125, "mm") ) +
  annotate(geom = "label", x = -4.85, y = 0.10, label = "KO+L", 
           colour = "white",  fill = "mediumpurple3", size = 2, 
           label.padding = grid::unit(0.5, "mm"), 
           label.r = unit(0.25, "mm"), family = "Arial",
           label.size = grid::unit(0.125, "mm") ) +
  # points on the right side
  geom_label_repel(data = subset(res_KOrS, p.adj < 0.1 & ratio > 0 ), 
                   aes(label = PTM), fill = 'white',
                   seed = 16, show.legend = F, segment.curvature = -1e-20, 
                   family = "Arial", size = 1.85, nudge_x = -0.25,
                   segment.color = 'black',verbose = F, 
                   box.padding = grid::unit(1, "mm"),
                   point.padding = grid::unit(0.55, "mm"),
                   label.padding = grid::unit(0.5, "mm"), 
                   label.size = 0.120, 
                   max.overlaps = 20)  +
  # points on the left side
  geom_label_repel(data = subset(res_KOrS,  p.adj < 0.1 & ratio < 0 ), 
                 aes(label = PTM), fill = 'white',
                 seed = 16, show.legend = F, segment.curvature = -1e-20, 
                 family = "Arial", size = 1.85, nudge_x = 0.25,
                 segment.color = 'black',verbose = F, 
                 box.padding = grid::unit(1, "mm"),
                 point.padding = grid::unit(0.55, "mm"),
                 label.padding = grid::unit(0.5, "mm"), 
                 label.size = 0.120, 
                 max.overlaps = 20)  +
  labs(x = expression(log[2] ~ "Fold Change"), y = expression(-log[10] ~ "P-Value")) +
  scale_fill_manual(values = c('dodgerblue', "gray84", "firebrick3"), guide = 'none') +
  guides(size = guide_legend(override.aes = list(fill = "white"))) +
  scale_x_continuous(expand = expansion(add = 0.02, mult = 0), limits = c(-5.5, 8.25), n.breaks = 7) +
  scale_y_continuous(expand = expansion(add = c(0, 0.05), mult = 0), limits = c(0, NA)) +
  scale_size_area(breaks = c(1e8, 1e9, 1e10, 1e11, 1e12),
                  labels = scientific_10, name = "WT PTM\nabundance") +
  theme_classic(base_size = 6, base_family = "Arial") +
  theme(panel.background = element_blank(),
        panel.border = element_blank(),
        panel.grid.major = element_line(linewidth = 0.2),
        legend.position = c(0.10, 0.83), 
        legend.key = element_blank(),
        legend.key.size = unit(1.0, units = 'mm'),
        legend.text = element_text(size = 5, margin = margin(l = -1, r = -1, unit = 'mm')),
        legend.direction = "vertical",
        legend.background = element_blank(),
        legend.box.background = element_blank(),
        legend.box.spacing = unit(-3, 'mm'),
        axis.text = element_text(colour = 'black'),
        axis.ticks.length = unit(1, units = 'mm'),
        axis.ticks = element_line(colour = 'black'),
        axis.title = element_text(size = 5),
        axis.line = element_line(linewidth = 0.2),
        plot.background = element_blank() ) -> p_rescues

p_rescues

Save to pdf.

Code
ggsave(path = pdf_dir_fig4, filename = "FigExtra_KOrS_vs_KOrL_Volcano_w_labels.pdf",
       plot = p_rescues, device = cairo_pdf, units = "cm",
       height = 6, width = 7.9)

I could still add:

  • rescue efficiency

  • heatmap with all PTMs.

Now plot to the supplementary figures

5 Supplementary Figure Panels

5.1 K27-K40 peptide modification ratios

5.1.1 H3.3K27

There’s no H3.3 signal in the main figure, here I plot only H3.3 signal.

Code
data_K27_36 |>
  group_by(Condition, K27, Histone) |>
  mutate(Average_Ratio = mean(Ratio)) |>
  summarise(Ratio_Sum = sum(Average_Ratio) / 3) |>
  subset(Histone %in% "H3.3") |>
  mutate(K27 = factor(K27, levels = c("K27ac", "K27me3", "K27me2", "K27me1", "unmod") ) ) |>
  ggplot() +
  aes(x = Condition, stratum = K27, alluvium = K27, 
      y = Ratio_Sum, fill = K27) +
  geom_alluvium(alpha = 0.5, width = 0.4) +
  geom_stratum2(alpha = 1, width = 0.6, colour = "black") +
  geom_fit_text(aes(label = K27 ),
                stat = "stratum", width = 0.6, min.size = 2,
                padding.x = grid::unit(0.2, "mm"),
                padding.y = grid::unit(0.2, "mm"),
                size = 7, family = "Arial", show.legend = F) +
  scale_fill_manual(values = c("K27ac" = "#feac81", "K27me1" ="#ced1af", 
                               "K27me2" = "#748f46", "K27me3" = "#47632a", 
                               "unmod" = "gray" ),
                    name = 'H3.3') +
  scale_x_discrete(expand = expansion(mult = 0.05, add = c(0, 0.1)),
                   labels = c('WT', 'CSex4', '∆ex4', 'KO', 
                              'KO+L', 'KO+S')) +
  scale_y_continuous(expand = expansion(mult = 0, add = 0) ) +
  guides(fill = guide_legend(byrow = TRUE) ) + # required for the legend item spacing 
  labs(y = "Average Ratio") +
  alluvial_thm -> p_Alluvial_H33K27

p_Alluvial_H33K27 

Save to pdf figure S4D.

Code
ggsave(filename = "FigS4D_H3.3_only_K27_Mrgd_Alluvial.pdf", path = pdf_dir_fig4, 
       plot = p_Alluvial_H33K27, device = cairo_pdf, units = "cm",
       width = 6.5, height = 4.85)

5.1.2 H3 and H3.3 K36

Plot the pileup of H3 & H3.3 K36 methylation.

Code
data_K27_36 |>
  group_by(Condition, K36, Histone) |>
  mutate(Average_Ratio = mean(Ratio)) |>
  summarise(Ratio_Sum = sum(Average_Ratio) / 3) |>
  subset(Histone == c("H3", "H3.3") ) |>
  mutate(K36 = factor(K36, levels = c("K36me3", "K36me2", "K36me1", "unmod") ) ) |>
  ggplot() +
  aes(x = Condition, stratum = K36, alluvium = K36, 
      y = Ratio_Sum, fill = K36) +
  facet_wrap(~ Histone, scales = 'free_y') +
  geom_alluvium(alpha = 0.5, width = 0.4) +
  geom_stratum(alpha = 1, width = 0.6, linewidth = 0.1, colour = "black") +
  geom_fit_text(aes(label = K36 ), stat = "stratum", width = 0.6, min.size = 2,
                padding.x = grid::unit(0.2, "mm"), size = 7, family = "Arial",
                padding.y = grid::unit(0.2, "mm"),
                show.legend = F) +
  scale_fill_manual(values = c("K36me1" = "#87B7AA", "K36me2" = "#8FA9C2",
                               "K36me3" = "#8386B7", "unmod" = "gray" ) ) +
  scale_x_discrete(expand = expansion(mult = 0.05, add = c(0, 0.1)),
                   labels = c('WT', 'CSex4', '∆ex4', 'KO', 'KO+L', 'KO+S')
                   ) +
  scale_y_continuous(expand = expansion(mult = 0, add = 0) ) +
  guides(fill = guide_legend(byrow = TRUE) ) + # required for the legend item spacing 
  labs(y = "Average Ratio") +
  alluvial_thm +
  theme(legend.title = element_blank()) -> p_Alluvial_K36

p_Alluvial_K36

Save to pdf figure S4E

Code
ggsave(path = pdf_dir_fig4, filename = "FigS4E_H3_H3.3_K36_Mrgd_Alluvial.pdf",
       plot = p_Alluvial_K36, device = cairo_pdf, units = "cm",
       width = 12.0, height = 5)

5.2 Extra plots

Here I plot extra histone PTMs that were not included in the publication in case someone is interested.

5.2.1 H3K4

Plot H3K4 modifications.

Code
data_dia_auto |>
  group_by(Condition, Modification) |>
  subset(Peptide_Start >= 3 & Peptide_End <= 8) |>
  mutate(Average_Ratio = mean(Ratio))  |>
  summarise(Ratio_Sum = sum(Average_Ratio) / 3, .groups = 'keep') |>
  ggplot() +
  aes(x = Condition, stratum = Modification, alluvium = Modification, 
      y = Ratio_Sum, fill = Modification) +
  geom_alluvium(alpha = 0.4, width = 0.4) +
  geom_stratum2(alpha = 1, width = 0.6, size = 0.1) +
  geom_fit_text(aes(label = Modification ),
                stat = "stratum", width = 0.6, min.size = 2,
                padding.x = grid::unit(0.2, "mm"),
                padding.y = grid::unit(0.2, "mm"),
                size = 7, family = "Arial", show.legend = F) +
  scale_fill_manual(values = met.brewer(name = "Demuth", direction = 1, n = 5), name = 'H3' ) +
  scale_x_discrete(expand = expansion(mult = 0.05, add = c(0, 0.1)),
                   labels = c('WT', 'CSex4', '∆ex4', 'KO', 
                              'KO+L', 'KO+S')) +
  scale_y_continuous(expand = expansion(mult = 0, add = 0) ) +
  guides(fill = guide_legend(byrow = TRUE) ) + # required for the legend item spacing 
  labs(y = "Average Ratio") +
  alluvial_thm +
  theme(legend.text = element_text(margin = margin(l = -9, unit = "mm")), 
        legend.key.width = unit(6.5, "mm") ) -> p_Alluvial_K4

p_Alluvial_K4

Save to pdf.

Code
ggsave(path = pdf_dir_fig4, filename = "FigExtra_H3K4_Mrgd_Alluvial.pdf",
       plot = p_Alluvial_K4, device = cairo_pdf, units = "cm",
       width = 6.5, height = 4.85)

5.2.2 H3K9

Plot all combinatorial K9 and K14 modifications on histone H3.

Code
data_dia_auto |>
  group_by(Condition, Modification) |>
  subset(Histone == "H3") |>
  subset(Peptide_Start >= 9 & Peptide_End <= 17) |>
  mutate(Average_Ratio = mean(Ratio) ) |>
  summarise(Ratio_Sum = sum(Average_Ratio) / 3, .groups = 'keep') |>
  ggplot() +
  aes(x = Condition, stratum = Modification, alluvium = Modification, 
      y = Ratio_Sum, fill = Modification) +
  geom_alluvium(alpha = 0.4, width = 0.4) +
  geom_stratum2(alpha = 1, width = 0.6, size = 0.1) +
  geom_fit_text(aes(label = Modification ),
                stat = "stratum", width = 0.6, min.size = 2,
                padding.x = grid::unit(0.2, "mm"),
                padding.y = grid::unit(0.2, "mm"),
                size = 7, family = "Arial", show.legend = F) +
  scale_fill_manual(values = met.brewer(name = "Hokusai1", direction = -1, n = 15), name = 'H3' ) + 
  scale_x_discrete(expand = expansion(mult = 0.05, add = c(0, 0.1)),
                   labels = c('WT', 'CSex4', '∆ex4', 'KO', 
                              'KO+L', 'KO+S')) +
  scale_y_continuous(expand = expansion(mult = 0, add = 0) ) +
  labs(y = "Average Ratio") +
  alluvial_thm +
  theme(legend.text = element_text(margin = margin(l = -12.0, unit = "mm")), 
        legend.key.width = unit(11.2, "mm"),
        legend.key.height = unit(1, 'mm'),
        legend.spacing.x = unit(1, 'mm') ) -> p_Alluvial_K9

p_Alluvial_K9

Code
ggsave(path = pdf_dir_fig4, filename = "FigExtra_H3K9_Mrgd_Alluvial.pdf",
       plot = p_Alluvial_K9, device = cairo_pdf, units = "cm",
       width = 6.4, height = 4.85)

5.2.3 All combinations of H3/H3.3 K27 & K36

Plot all combinatorial K27 and K36 modifications on histone H3 and H3.3 variant.

Code
data_dia_auto |>
  group_by(Condition, Modification, Histone) |>
  subset(Peptide_Start >= 27 & Peptide_End <= 40) |>
  mutate(Average_Ratio = mean(Ratio))  |>
  summarise(Ratio_Sum = sum(Average_Ratio) / 3, .groups = 'keep') |>
  ggplot() +
  aes(x = Condition, stratum = Modification, alluvium = Modification, 
      y = Ratio_Sum, fill = Modification) +
  facet_wrap(~ Histone) +
  geom_alluvium(alpha = 0.4, width = 0.4) +
  geom_stratum(alpha = 1, width = 0.6, size = 0.1) +

  geom_fit_text(aes(label = Modification ),
                stat = "stratum", width = 0.6, min.size = 2,
                padding.x = grid::unit(0.2, "mm"),
                padding.y = grid::unit(0.2, "mm"),
                size = 7, family = "Arial", show.legend = F) +
  scale_fill_manual(values = met.brewer(name = "Hiroshige", direction = 1, n = 15) ) +
  scale_x_discrete(expand = expansion(mult = 0.05, add = c(0, 0.1)),
                   labels = c('WT', 'CSex4', '∆ex4', 'KO', 
                              'KO+L', 'KO+S'), name = 'H3 & H3.3') +
  scale_y_continuous(expand = expansion(mult = 0, add = 0) ) +
  labs(y = "Average Ratio") +
  alluvial_thm +
  theme(legend.text = element_text(margin = margin(l = -14.5, unit = "mm")), 
        legend.key.width = unit(13.75, "mm"),
        legend.key.height = unit(1, 'mm'),
        legend.spacing.x = unit(1, 'mm') ) -> p_Alluvial_K27_K36

p_Alluvial_K27_K36

This is the best representative picture of K27-K40 modifications relative abundance. This also confirms how the H3.3K36me3 is really the most dramatic change between KO+L and KO+S.

Code
ggsave(path = pdf_dir_fig4, filename = "FigExtra_H3_H3.3_K27_K36_Mrgd_Alluvial.pdf",
       plot = p_Alluvial_K27_K36, device = cairo_pdf, units = "cm",
       width = 14, height = 5)

End analysis.

5.3 Protein PSI quantification

Code
Protein_PSI_3reps_LS <- read_excel(path = protein_psi_path, sheet = "Sheet1") |>
  pivot_longer(cols = starts_with('SUZ12'), names_to = 'Proteoform', values_to = 'PSI') |>
  mutate(PSI = PSI * 100)

Protein_PSI_3reps_LS |> group_by(Proteoform) |> summarise(median(PSI))
# A tibble: 2 × 2
  Proteoform `median(PSI)`
  <chr>              <dbl>
1 SUZ12-L            90.6 
2 SUZ12-S             9.43

Plot

Code
ggplot(Protein_PSI_3reps_LS) +
  aes(x = Proteoform, y = PSI, fill = Proteoform) +
  geom_boxplot(show.legend = F, outlier.shape = NA, width = 0.5, lwd = 0.1) +
  geom_point(shape = 21, size = 0.85, stroke = 0.1, show.legend = F) +
  coord_cartesian(ylim = c(0, 100), clip = 'off') +
  scale_fill_manual(values = c("SUZ12-L" = "mediumpurple3", "SUZ12-S" = "darkorange1")) + 
  scale_y_continuous(expand = expansion(mult = 0, add = 0),
                     n.breaks = 10) +
  labs(x = 'Proteoform', 
       y = 'Rel. abundance in WT ESCs') +
  alluvial_thm +
  theme(panel.grid.major.y = element_line(colour = 'gray84', linewidth = 0.2),
        axis.title.x = element_text(size = 5) ) -> p_prot_PSI
p_prot_PSI

Save to pdf figure S4C

Code
ggsave(filename = "FigS4C_SUZ12_WT_proteoforms_quantification_barplot.pdf", 
       plot = p_prot_PSI, device = cairo_pdf, path = pdf_dir_fig4, units = 'cm',
       width = 2.65, height = 2.75)

6 Session Info

Code
sessioninfo::session_info()
─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.1.2 (2021-11-01)
 os       macOS Catalina 10.15.7
 system   x86_64, darwin17.0
 ui       X11
 language (EN)
 collate  C
 ctype    UTF-8
 tz       Europe/Madrid
 date     2024-02-07
 pandoc   2.10.1 @ /usr/local/bin/ (via rmarkdown)

─ Packages ───────────────────────────────────────────────────────────────────
 ! package              * version   date (UTC) lib source
   ade4                   1.7-22    2023-02-06 [1] CRAN (R 4.1.2)
   affy                   1.72.0    2021-10-26 [1] Bioconductor
   affyio                 1.64.0    2021-10-26 [1] Bioconductor
   annotate               1.72.0    2021-10-26 [1] Bioconductor
   AnnotationDbi          1.56.2    2021-11-09 [1] Bioconductor
   assertthat             0.2.1     2019-03-21 [1] CRAN (R 4.1.0)
   Biobase                2.54.0    2021-10-26 [1] Bioconductor
   BiocFileCache          2.2.1     2022-01-23 [1] Bioconductor
   BiocGenerics           0.40.0    2021-10-26 [1] Bioconductor
   BiocManager            1.30.22   2023-08-08 [1] CRAN (R 4.1.2)
   BiocParallel           1.28.3    2021-12-09 [1] Bioconductor
   biomaRt                2.50.3    2022-02-06 [1] Bioconductor
   Biostrings             2.62.0    2021-10-26 [1] Bioconductor
   bit                    4.0.5     2022-11-15 [1] CRAN (R 4.1.2)
   bit64                  4.0.5     2020-08-30 [1] CRAN (R 4.1.0)
   bitops                 1.0-7     2021-04-24 [1] CRAN (R 4.1.0)
   blob                   1.2.4     2023-03-17 [1] CRAN (R 4.1.2)
   bslib                  0.5.1     2023-08-11 [1] CRAN (R 4.1.2)
   cachem                 1.0.8     2023-05-01 [1] CRAN (R 4.1.2)
   Cairo                * 1.6-1     2023-08-18 [1] CRAN (R 4.1.2)
   callr                  3.7.3     2022-11-02 [1] CRAN (R 4.1.2)
   cellranger             1.1.0     2016-07-27 [1] CRAN (R 4.1.0)
   circlize               0.4.15    2022-05-10 [1] CRAN (R 4.1.2)
   cli                    3.6.1     2023-03-23 [1] CRAN (R 4.1.2)
   clue                   0.3-64    2023-01-31 [1] CRAN (R 4.1.2)
   cluster                2.1.4     2022-08-22 [1] CRAN (R 4.1.2)
   codetools              0.2-19    2023-02-01 [1] CRAN (R 4.1.2)
   colorspace             2.1-0     2023-01-23 [1] CRAN (R 4.1.2)
   ComplexHeatmap         2.10.0    2021-10-26 [1] Bioconductor
   crayon                 1.5.2     2022-09-29 [1] CRAN (R 4.1.2)
   crosstalk              1.2.0     2021-11-04 [1] CRAN (R 4.1.0)
   csaw                   1.28.0    2021-10-26 [1] Bioconductor
   curl                   5.2.0     2023-12-08 [1] CRAN (R 4.1.2)
   DBI                    1.1.3     2022-06-18 [1] CRAN (R 4.1.2)
   dbplyr                 2.3.3     2023-07-07 [1] CRAN (R 4.1.2)
   DelayedArray           0.20.0    2021-10-26 [1] Bioconductor
   DEP                  * 1.16.0    2021-10-26 [1] Bioconductor
   desc                   1.4.2     2022-09-08 [1] CRAN (R 4.1.2)
   DESeq2                 1.34.0    2021-10-26 [1] Bioconductor
   devtools               2.4.5     2022-10-11 [1] CRAN (R 4.1.2)
   digest                 0.6.33    2023-07-07 [1] CRAN (R 4.1.2)
   doParallel             1.0.17    2022-02-07 [1] CRAN (R 4.1.2)
   dplyr                * 1.1.2     2023-04-20 [1] CRAN (R 4.1.2)
   DT                   * 0.29      2023-08-29 [1] CRAN (R 4.1.2)
   edgeR                  3.36.0    2021-10-26 [1] Bioconductor
   ellipsis               0.3.2     2021-04-29 [1] CRAN (R 4.1.0)
   evaluate               0.21      2023-05-05 [1] CRAN (R 4.1.2)
   fansi                  1.0.4     2023-01-22 [1] CRAN (R 4.1.2)
   farver                 2.1.1     2022-07-06 [1] CRAN (R 4.1.2)
   fastmap                1.1.1     2023-02-24 [1] CRAN (R 4.1.2)
   fdrtool                1.2.17    2021-11-13 [1] CRAN (R 4.1.0)
   filelock               1.0.2     2018-10-05 [1] CRAN (R 4.1.0)
   forcats                1.0.0     2023-01-29 [1] CRAN (R 4.1.2)
   foreach                1.5.2     2022-02-02 [1] CRAN (R 4.1.2)
   foreign                0.8-84    2022-12-06 [1] CRAN (R 4.1.2)
   fs                     1.6.3     2023-07-20 [1] CRAN (R 4.1.2)
   genefilter             1.76.0    2021-10-26 [1] Bioconductor
   geneplotter            1.72.0    2021-10-26 [1] Bioconductor
   generics               0.1.3     2022-07-05 [1] CRAN (R 4.1.2)
   GenomeInfoDb           1.30.1    2022-01-30 [1] Bioconductor
   GenomeInfoDbData       1.2.7     2023-08-29 [1] Bioconductor
   GenomicRanges          1.46.1    2021-11-18 [1] Bioconductor
   GetoptLong             1.0.5     2020-12-15 [1] CRAN (R 4.1.0)
   ggalluvial             0.12.5    2023-02-22 [1] CRAN (R 4.1.2)
   ggfittext              0.10.0    2023-04-04 [1] CRAN (R 4.1.2)
   ggplot2                3.4.3     2023-08-14 [1] CRAN (R 4.1.2)
   ggrepel                0.9.3     2023-02-03 [1] CRAN (R 4.1.2)
   ggseqlogo              0.1       2017-07-25 [1] CRAN (R 4.1.0)
   GlobalOptions          0.1.2     2020-06-10 [1] CRAN (R 4.1.0)
   glue                   1.6.2     2022-02-24 [1] CRAN (R 4.1.2)
   gmm                    1.8       2023-06-06 [1] CRAN (R 4.1.2)
   gtable                 0.3.4     2023-08-21 [1] CRAN (R 4.1.2)
   hexbin                 1.28.3    2023-03-21 [1] CRAN (R 4.1.2)
   hms                    1.1.3     2023-03-21 [1] CRAN (R 4.1.2)
   htmltools              0.5.6     2023-08-10 [1] CRAN (R 4.1.2)
   htmlwidgets            1.6.2     2023-03-17 [1] CRAN (R 4.1.2)
   httpuv                 1.6.11    2023-05-11 [1] CRAN (R 4.1.2)
   httr                   1.4.7     2023-08-15 [1] CRAN (R 4.1.2)
   impute                 1.68.0    2021-10-26 [1] Bioconductor
   imputeLCMD             2.1       2022-06-10 [1] CRAN (R 4.1.2)
   IRanges                2.28.0    2021-10-26 [1] Bioconductor
   iterators              1.0.14    2022-02-05 [1] CRAN (R 4.1.2)
   jquerylib              0.1.4     2021-04-26 [1] CRAN (R 4.1.0)
   jsonlite               1.8.7     2023-06-29 [1] CRAN (R 4.1.2)
   KEGGREST               1.34.0    2021-10-26 [1] Bioconductor
   knitr                  1.43      2023-05-25 [1] CRAN (R 4.1.2)
   labeling               0.4.2     2020-10-20 [1] CRAN (R 4.1.0)
   later                  1.3.1     2023-05-02 [1] CRAN (R 4.1.2)
   lattice                0.21-8    2023-04-05 [1] CRAN (R 4.1.2)
   lifecycle              1.0.3     2022-10-07 [1] CRAN (R 4.1.2)
   limma                  3.50.3    2022-04-07 [1] Bioconductor
   locfit                 1.5-9.8   2023-06-11 [1] CRAN (R 4.1.2)
   magick                 2.7.5     2023-08-07 [1] CRAN (R 4.1.2)
   magrittr               2.0.3     2022-03-30 [1] CRAN (R 4.1.2)
   MALDIquant             1.22.1    2023-03-20 [1] CRAN (R 4.1.2)
   MASS                   7.3-60    2023-05-04 [1] CRAN (R 4.1.2)
   Matrix                 1.6-1     2023-08-14 [1] CRAN (R 4.1.2)
   MatrixGenerics         1.6.0     2021-10-26 [1] Bioconductor
   matrixStats            1.0.0     2023-06-02 [1] CRAN (R 4.1.2)
   memoise                2.0.1     2021-11-26 [1] CRAN (R 4.1.0)
   metapod                1.2.0     2021-10-26 [1] Bioconductor
   MetBrewer              0.2.0     2022-03-21 [1] CRAN (R 4.1.2)
   mime                   0.12      2021-09-28 [1] CRAN (R 4.1.0)
   miniUI                 0.1.1.1   2018-05-18 [1] CRAN (R 4.1.0)
   msa                    1.26.0    2021-10-26 [1] Bioconductor
   MsCoreUtils            1.6.2     2022-02-24 [1] Bioconductor
   MSnbase                2.20.4    2022-01-16 [1] Bioconductor
   munsell                0.5.0     2018-06-12 [1] CRAN (R 4.1.0)
   mvtnorm                1.2-3     2023-08-25 [1] CRAN (R 4.1.2)
   mzID                   1.32.0    2021-10-26 [1] Bioconductor
   mzR                    2.28.0    2021-10-27 [1] Bioconductor
   ncdf4                  1.21      2023-01-07 [1] CRAN (R 4.1.2)
 R niar                 * 0.3.0     <NA>       [?] <NA>
   norm                   1.0-11.1  2023-06-18 [1] CRAN (R 4.1.2)
   patchwork              1.1.3     2023-08-14 [1] CRAN (R 4.1.2)
   pcaMethods             1.86.0    2021-10-26 [1] Bioconductor
   pheatmap             * 1.0.12    2019-01-04 [1] CRAN (R 4.1.0)
   pillar                 1.9.0     2023-03-22 [1] CRAN (R 4.1.2)
   pkgbuild               1.4.2     2023-06-26 [1] CRAN (R 4.1.2)
   pkgconfig              2.0.3     2019-09-22 [1] CRAN (R 4.1.0)
   pkgload                1.3.2.1   2023-07-08 [1] CRAN (R 4.1.2)
   plyr                   1.8.8     2022-11-11 [1] CRAN (R 4.1.2)
   png                    0.1-8     2022-11-29 [1] CRAN (R 4.1.2)
   preprocessCore         1.56.0    2021-10-26 [1] Bioconductor
   prettyunits            1.1.1     2020-01-24 [1] CRAN (R 4.1.0)
   processx               3.8.2     2023-06-30 [1] CRAN (R 4.1.2)
   profvis                0.3.8     2023-05-02 [1] CRAN (R 4.1.2)
   progress               1.2.2     2019-05-16 [1] CRAN (R 4.1.0)
   promises               1.2.1     2023-08-10 [1] CRAN (R 4.1.2)
   ProtGenerics           1.26.0    2021-10-26 [1] Bioconductor
   ps                     1.7.5     2023-04-18 [1] CRAN (R 4.1.2)
   psychTools             2.3.6     2023-06-18 [1] CRAN (R 4.1.2)
   purrr                  1.0.2     2023-08-10 [1] CRAN (R 4.1.2)
   R6                     2.5.1     2021-08-19 [1] CRAN (R 4.1.0)
   rappdirs               0.3.3     2021-01-31 [1] CRAN (R 4.1.0)
   RColorBrewer           1.1-3     2022-04-03 [1] CRAN (R 4.1.2)
   Rcpp                   1.0.11    2023-07-06 [1] CRAN (R 4.1.2)
   RCurl                  1.98-1.12 2023-03-27 [1] CRAN (R 4.1.2)
   readr                * 2.1.4     2023-02-10 [1] CRAN (R 4.1.2)
   readxl               * 1.4.3     2023-07-06 [1] CRAN (R 4.1.2)
   remotes                2.4.2.1   2023-07-18 [1] CRAN (R 4.1.2)
   rjson                  0.2.21    2022-01-09 [1] CRAN (R 4.1.2)
   rlang                  1.1.1     2023-04-28 [1] CRAN (R 4.1.2)
   rmarkdown              2.24      2023-08-14 [1] CRAN (R 4.1.2)
   rprojroot              2.0.3     2022-04-02 [1] CRAN (R 4.1.2)
   Rsamtools              2.10.0    2021-10-26 [1] Bioconductor
   RSQLite                2.3.1     2023-04-03 [1] CRAN (R 4.1.2)
   rstudioapi             0.15.0    2023-07-07 [1] CRAN (R 4.1.2)
   S4Vectors              0.32.4    2022-03-29 [1] Bioconductor
   sandwich               3.0-2     2022-06-15 [1] CRAN (R 4.1.2)
   sass                   0.4.7     2023-07-15 [1] CRAN (R 4.1.2)
   scales                 1.2.1     2022-08-20 [1] CRAN (R 4.1.2)
   seqinr                 4.2-30    2023-04-05 [1] CRAN (R 4.1.2)
   sessioninfo            1.2.2     2021-12-06 [1] CRAN (R 4.1.0)
   shape                  1.4.6     2021-05-19 [1] CRAN (R 4.1.0)
   shiny                  1.7.5     2023-08-12 [1] CRAN (R 4.1.2)
   shinydashboard         0.7.2     2021-09-30 [1] CRAN (R 4.1.0)
   stringi                1.7.12    2023-01-11 [1] CRAN (R 4.1.2)
   stringr                1.5.0     2022-12-02 [1] CRAN (R 4.1.2)
   SummarizedExperiment   1.24.0    2021-10-26 [1] Bioconductor
   survival               3.5-7     2023-08-14 [1] CRAN (R 4.1.2)
   tibble                 3.2.1     2023-03-20 [1] CRAN (R 4.1.2)
   tidyr                  1.3.0     2023-01-24 [1] CRAN (R 4.1.2)
   tidyselect             1.2.0     2022-10-10 [1] CRAN (R 4.1.2)
   tmvtnorm               1.5       2022-03-22 [1] CRAN (R 4.1.2)
   tzdb                   0.4.0     2023-05-12 [1] CRAN (R 4.1.2)
   urlchecker             1.0.1     2021-11-30 [1] CRAN (R 4.1.0)
   usethis                2.2.2     2023-07-06 [1] CRAN (R 4.1.2)
   utf8                   1.2.3     2023-01-31 [1] CRAN (R 4.1.2)
   vctrs                  0.6.3     2023-06-14 [1] CRAN (R 4.1.2)
   vroom                  1.6.3     2023-04-28 [1] CRAN (R 4.1.2)
   vsn                    3.62.0    2021-10-26 [1] Bioconductor
   withr                  2.5.0     2022-03-03 [1] CRAN (R 4.1.2)
   xfun                   0.40      2023-08-09 [1] CRAN (R 4.1.2)
   XICOR                  0.4.1     2023-04-21 [1] CRAN (R 4.1.2)
   XML                    3.99-0.14 2023-03-19 [1] CRAN (R 4.1.2)
   xml2                   1.3.5     2023-07-06 [1] CRAN (R 4.1.2)
   xtable                 1.8-4     2019-04-21 [1] CRAN (R 4.1.0)
   XVector                0.34.0    2021-10-26 [1] Bioconductor
   yaml                   2.3.7     2023-01-23 [1] CRAN (R 4.1.2)
   zlibbioc               1.40.0    2021-10-26 [1] Bioconductor
   zoo                    1.8-12    2023-04-13 [1] CRAN (R 4.1.2)

 [1] /Library/Frameworks/R.framework/Versions/4.1/Resources/library

 R ── Package was removed from disk.

──────────────────────────────────────────────────────────────────────────────

References

Yuan, Zuo-Fei, Simone Sidoli, Dylan M Marchione, Johayra Simithy, Kevin A Janssen, Mary R Szurgot, and Benjamin A Garcia. 2018. EpiProfile 2.0: A Computational Platform for Processing Epi-Proteomics Mass Spectrometry Data.” Journal of Proteome Research 17 (7): 2533–41. https://doi.org/10.1021/acs.jproteome.8b00133.
Zhang, Xiaofei, Arne H. Smits, Gabrielle BA van Tilburg, Huib Ovaa, Wolfgang Huber, and Michiel Vermeulen. 2018. Proteome-wide identification of ubiquitin interactions using UbIA-MS.” Nature Protocols 13 (3): 530–50. https://doi.org/10.1038/nprot.2017.147.